Algorithm for molecular dynamics simulations of spin liquids 
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A new symplectic time-reversible algorithm for numerical integration of the equations of motion 
in magnetic liquids is proposed. It is tested and applied to molecular dynamics simulations of 
a Heisenberg spin fluid. We show that the algorithm exactly conserves spin lengths and can be 
used with much larger time steps than those inherent in standard predictor-corrector schemes. The 
results obtained for time correlation functions demonstrate the evident dynamic interplay between 
the liquid and magnetic subsystems. 
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Computer experiments remain an important tool for 
the prediction and theoretical understanding of vari- 
ous phenomena in magnetic materials. The methods of 
Monte-Carlo (MC) and molecular dynamics (MD) were 
intensively exploited over the years for the investigation 
of phase diagrams, critical phenomena, scaling, and the 
dynamic behavior of lattice systems such as the Ising, the 
XY, and the Heisenberg model [jl] [| . 

The necessity to extend these studies to disordered 
models of magnetic liquids was motivated by a great 
amount of additional physical properties arising when 
both spin (orientational) and liquid (translational) de- 
grees of freedom are taken into account |4| |]]. The com- 
puter experiments for such systems have been restricted 
to MC simulations [j|,f7| in which only static quantities 
could be calculated. Dynamic phenomena, in particular, 
spin and density relaxations, and the effects connected 
with the mutual influence of magnetic and liquid subsys- 
tems can be investigated in MD simulations. Our special 
interest to this problem was also stimulated by the re- 
sults obtained recently for a Heisenberg fluid within the 
hydrodynamic theory |j] . One of them was the predic- 
tion that the shape of magnetic dynamic structure factor 
can change qualitatively in comparison with the lattice 
model due to the coupling between the subsystems. 

Until now, there have been no attempts to simulate 
spin liquids within the MD approach. This can be ex- 
plained by the absence of an efficient MD algorithm for 
handling the corresponding equations of motion (EOM). 
The traditional numerical methods |10| for solving dif- 
ferential equations are unsuitable because they become 
highly unstable on time scales used in MD simulations. 
As has been well established for pure liquid systems 
jll],[L2]], even standard predictor-corrector schemes are 
not efficient because of poor total energy conservation. 

The properties of an acceptable algorithm for long- 
time observations over a many-body system should be: 
stability, accuracy, speed and ease of implementation. 
There exists only a small group of integrators satisfying 
these criteria. An important one is the velocity Verlet 
(VV) algorithm (l^Jl^l which allows a high accuracy with 



minimal costs in terms of time-consuming force evalua- 
tions. However, the VV and other similar schemes [|ll|jl5| 
were designed to simulate pure liquid dynamics. In the 
case of magnetic liquids the situation is more complicated 
since the translational positions and momenta are coupled 
with spin orientations in a characteristic way and, hence, 
all these dynamical variables must be considered simul- 
taneously. This requires substantial revision of the liquid 
dynamic algorithms. 

Recently, new algorithms have been devised for spin 
dynamics simulations of lattice systems [ fl6| . They are 
based (like the VV integrator) on the Suzuki- Trotter (ST) 
decomposition method and appear to be much more effi- 
cient than predictor-corrector schemes. These algorithms 
are applicable to spin systems if the decomposition on 
two (or several) noninteracting sublattices is possible. 
However, they cannot be used for models with arbitrary 
spatial spin distributions and, therefore, not for spin liq- 
uids. 

In this Letter we develop the idea of using ST-like de- 
compositions for spin liquid dynamics and derive the de- 
sired MD algorithm. This allows quantitative measure- 
ment of dynamical structure factors of a Heisenberg fluid. 
The main result obtained (reflecting the influence of the 
liquid sybsystem on spin dynamics) is the identification 
of a new propagative sound-like mode in the spectrum of 
collective longitudinal spin excitations. 

Consider a classical system composed of N magnetic 
particles of mass m, described by the Hamiltonian |j],U 
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Here and Vj are the translational position and veloc- 
ity, respectively, of particle i carrying spin s^. The liquid 
potential is denoted by $(ry), and J{rij) > is the 
exchange integral for a pair of spins with interparticle 
distance r,j . The classical approach treats Sj as a three- 
component continuous vector with a fixed length for each 
site i. We put for convenience |sj| = 1, so that J is mea- 
sured in energy units. 
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In order to study the dynamic properties, the equations 
of motion given by dp/dt = Lp(t) must be integrated nu- 
merically, where 
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= ^2( L n + L v, +L Si ) =L r + L v + L s , (2) 
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is the Liouville operator, i.e., Lp = [p,H] with [, ] be- 
ing the Poisson bracket; p = {r.^v^Sj} denotes the 
full set of microscopic phase variables; = fj/m and 
u)i — —gi/h are the acceleration and local Larmor fre- 
quency, respectively, with = -J2j(j 7 ti)[d®( r ij)/drij - 
dJ(ry)/dry s l -s j ]v ij /r i: j and g, = E j(j#i) Jir./is, being 
the force and internal magnetic field. Note that the op- 
erators L r , L v and L s act only on position, velocity and 
spin, respectively, and the quantum Poisson bracket was 
used to obtain L s . 

The solutions can be cast in the form p{t + h) = 
e Lh p(t) = e ( - Lr+Lv+L ' )h p(t), where h is the time step. 
Since the exponential propagator e Lh cannot be evalu- 
ated exactly, one introduces some approximations which 
take advantage of the smallness of h. Assuming for 
the moment that spin variables are frozen, i.e. set- 
ting L s — > 0, we come to the usual (liquid-like) EOM. 
They can be solved in a quite efficient way using the 
second-order VV integrator [ |13|Jl4[| which is based on 
the ST formula e ^+ L ^ h = + Qfj^y 

Taking into account the fact that this formula is valid 
for arbitrary two operators and unfreezing now the 
spin variables, we obtain immediately: e ( ir+iv+Ls )' 1 = 
e L v h/2 e (L r +L e )h e L v h/2 + o(h 3 ), where the sum L r + L s 

was treated as one operator. The spin-position sub- 
propagator can further be decomposed in a similar way, 
c (L r +L B )h = e L r h/2 e L B h e L r h/2 + o(h 3 ), resulting in a full 

propagation of the form 



p(t + h) = e L ^e L ^e L * h e L ^e L ^p{t) + 0{h 3 ) . (3) 

Note that other decompositions are also possible, but 
then the local fields gi and/or forces fi have to be updated 
(the most time-consuming operations) more frequently, 
which reduces the efficiency of the computations. 

The main idea of the decompositions is to obtain sub- 
propagators which can be evaluated analytically. It can 
be shown H] that the position e LrT = Y[i eL " iT an d ve ~ 
locity e LvT = ]X e iv i T propagations represent shift oper- 
ators, namely, e r » r ri = r s ; + v^r and e Lv i T Vi = Vj +a;T. 
Since the components L r and L v (as well as L s ) do not 
commute, such shifts must be performed in a rigorous 
order (as specified by Eq. (3)) and applied to the current 
values of and v, within the time step. 

The spin subdynamics is described in Eq. (3) by the 
exponential operator e L " h . This operator has no sim- 
ple explicit form, because the Larmor frequency u>j for 



each particle depends in general on the orientations of 
all other spins of the system. The explicit solution, nev- 
ertheless, may be found as follows. Since all the partial 
components L Si do not commute each other, it is quite 
natural to find an ST-like decomposition for the whole 
set of these operators. This results to the expression 
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which constitutes an ST analog for arbitrary number of 
operators and is accurate to the same order 0(h 3 ) as the 
terms already truncated. Again, other 0(h 3 ) decomposi- 
tions may be introduced. However, only Eq. (4) will lead 
to a scheme with minimal number of local field recalcu- 
lations. 

The problem is now considerably simplified because, 
according to Eq. (4), each current value of s, is updated 
spin by spin at a fixed instantaneous Larmor frequency 
uji, and this case allows analytical solutions: e L *i T Si(t) = 
Bi(t,T)si(t). Here Di(t,r) = I + W, sin(^r) + Wf(l - 
cos(wir)) denotes an orthonormal (DD + = I) matrix 
of rotation around axis u)i on angle u>iT and W, = 
W(u>j) is a skewsymmetric matrix (W a p = — W^ a ) 

with *WxY = -W^, W^z = LUy, = an d 

u) = u)/lj. Since the decompositions used are correct 
within an uncertainty of order 0(h 3 ), the trigonometric 
functions can be replaced by their rational counterparts 
(see, e.g., @), cos£ = (1 - £ 2 /T)/(l + C 2 /4) + 0(£ 3 ) 
and sin £ = + ^ 2 / 4 ) + C(C 3 ), which maintain the 
orthonormality of D and are more efficient for the com- 
putations. Then the spin rotation reduces to 



* T Si(i) = \ si(t) + [wiXSi(i)]T+ —(u}i(u)i-Si(t)) 
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This completes the new algorithm. 

We note that our basic EOM are time-reversible and 
exact solutions behave symplectically. As can be shown, 
the algorithm derived reproduces these features, even 
though the trajectories are generated with a limited ac- 
curacy. Indeed, the initial propagator was decomposed 
(Eqs. (3) and (4)) into subparts symmetrically, and, as 
a consequence, the final expressions for rj, Vj and will 
be invariant with respect to the transformation h — > —h. 
Furthermore, simple shifts (applied separately in posi- 
tion and velocity space) do not change the phase volume. 
These properties are very important for our purpose be- 

. ,14|], the stability of 
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cause, as is now well established 

an algorithm normally follows from its time reversibility 
and symplecticity. Another nice property of the algo- 
rithm is its exact conservation of spin lengths (rotations 
given by Eq. (5) do not change the norm of vectors) that 
is crucial for the class of models considered. 

In our MD study of the Heisenberg fluid, we have used 
the Yukawa potential 0, J(r) = {ea/r) exp[(er — r)/er], 
and a soft-core potential pT| , $(r) = Ae[{a/r) 12 — 
(cr/r) 6 ] + e at r < 2 1 / 6 er and $(r) = otherwise, 
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for the description of spin and liquid interactions with 
the intensities e and e, respectively. The function J(r) 
was truncated at R = 2.5a and shifted to be zero at 
the truncation point to avoid force singularities. The 
simulations were carried out for N = 1000 particles 
(employing periodic boundary conditions) at a reduced 
density n* — Na 3 /V — 0.6, a reduced temperature 
T* = 'k B T/e = 1.5 < T* (where T* ~ 2.06 is the tem- 
perature of ferromagnetic transition a reduced core 
intensity e/e = 1, and a dynamical coupling parameter 
d = ff(me) 1 / 2 /'i = 2. This last parameter presents, 
in fact, the ratio r tr /r s , where r tr = o~{m/e) 1 / 2 and 
t s = h/e are the characteristic time intervals of vary- 
ing translational and spin variables, respectively. Since 
we are investigating a ferrophase and dealing with a mi- 
crocanonical (NVES) ensemble, a non-zero magnetiza- 
tion of the system must be specified additionally. This 
quantity was taken from our single MC simulation [fTs] , 
(S) /iV = 0.6536±0.0001, where ( ) denotes the canoni- 
cal averaging. All test runs were started from an identical 
well equilibrated configuration. The recalculation of lo- 
cal magnetic fields (during spin subdynamics (4)) took 
approximately the same processor time as that of trans- 
lational forces, spending in total 0.5 sec per step on the 
Origin 2000 workstation. It is worth emphasizing that 
contrary to spin lattice dynamics iful (when auxiliary 
MC cycles are involved to generate equilibrium configura- 
tions as initial conditions for the EOM), the equilibration 
of our system can be performed within NVES MD simu- 
lations exclusively (at the specified value for S = (S} ). 
This is possible because of the energy exchange between 
the spin and liquid subsystems. 

Symmetries of Hamiltonian (1) impose conservation 
laws on the total momentum P = rn^2 i \i, total spin 
S = J2i s i an d total energy E = H. These three inte- 
grals of motion cannot be conserved perfectly at the same 
time within any approximate scheme known. This a typ- 
ical situation in MD simulations. The MD results for the 
total energy E* = E/e (subsets (a)-(d)) and total spin S 
(subsets (e)-(h)) as functions of the length of the simu- 
lations are presented in Fig. 1. Four time steps, namely, 
h* = ft/r tr = 0.00125, 0.0025, 0.005 and 0.01, were used 
to integrate the EOM (solid curves). These results are 
compared with those obtained by us using the well es- 
tablished Adams-Bashforth-Moulton (ABM) predictor- 
corrector scheme |lCj (dashed curves in subsets (a) and 
(b)). As can be seen from Fig. la, the ABM integrator 
fulfills energy conservation up to a similar accuracy as our 
algorithn at the smallest time step h* = 0.00125. How- 
ever, for larger step sizes (see Fig. lb) the ABM scheme is 
unstable and, thus, cannot be used. Note that very small 
step sizes are impractical because then too much time- 
consuming force and field evaluations have to be done 
during the typical observation times. 

No systematic drift in E(t) and S(t) was observed 
within our algorithm at time steps up to h* = 0.01 
over a length of t/h = 100 000. The precision of the 
algorithm was measured in terms of the ratio Te — 



(((E(t) - E(0)) 2 )/({U{t) - [/(0)) 2 }) 1 / 2 of total and po- 
tential (U) energy fluctuations. Taking into account that 
for our system ((U(t) - CZ(O)) 2 ) 1 / 2 w 0.0335, we have ob- 
tained: T E w 0.12%, 0.28%, 0.98% and 7.7% for the 
time steps ft* = 0.00125, 0.0025, 0.005 and 0.01, respec- 
tively. In order to reproduce properly the features of mi- 
crocanonical ensembles the ratio Te should not exceed 
a few per cent. As we can see, time steps of ft* < 0.01 
satisfy this requirement and, thus, they can be used for 
precise calculations. 



FIG. 1. Reduced total energy E*(t)/N [(a)-(d)] and mag- 
netization S(t)/N [(e)-(h)] per spin as functions of the obser- 
vation time obtained within the decomposition [solid curves] 
and predictor-corrector [dashed curves in (a) and (b)] algo- 
rithms at four fixed time steps: h* = 0.00125, 0.0025, 0.005, 
and 0.01. The values E*(0)/N and 5(0) /N are plotted by the 
horizontal thin lines. 



Note that the decomposition and ABM methods con- 
serve the total momentum P to within machine accuracy. 
The reason is that all velocities are updated simultane- 
ously and the interparticle forces are evaluated exploiting 
Newton's third law. For similar reasons, the ABM inte- 
gration maintains the total magnetization S (but it does 
not conserve spin lengths). In our scheme the magnetiza- 
tion is not conserved exactly. However, the fluctuations 
appear to be very small (see Fig. le-lh) and lead to the 
values «(S(i) - S(O)) 2 }) 1 / 2 ~ 10" 7 , 5 • 10~ 7 , 2 • 10~ 6 and 
10" 5 at ft* =0.00125, 0.0025, 0.005 and 0.01, respectively. 

The spectra F(k,co) = ^ F(k, t)e~ iuJt dt of the 
spin-spin F^ T (k,t) = s^ T (0)- 8 J' T (t)ny(k, *)), 

density-density F nn (k,t) = . riij (k,t)) and spin- 

density Fl n {k : t) = <E M sH0K('M)> = F*u{k,t) time 
correlation functions are shown in Fig. 2. The super- 
scripts (L) and (T) refer to the longitudinal and trans- 
verse components of Sj with respect to the vector S, and 
nij(k,t) = N~ x exp{ik-(ri(0) — Tj(t))}. These functions 
were obtained within our decomposition integration at 
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h* = 0.005, and the microcanonical averaging ( ) was 
taken over 100 000 steps for each of 10 independent 
runs. A dimensionless representation has been used for 
F*(k,u>) = F(k,u>)/Tti with lu* — u7T tr , k* = fc//c m i n and 
fe min = 2^/^/3. 



In conclusion, we list the chief advantages of the new 
algorithm over existing numerical schemes: (i) time- 
reversibility and symplecticity; (ii) explicitness (no iter- 
ation); (iii) exact conservation of spin lengths; (iv) much 
more accuracy in total energy conservation. Moreover, 
its excellent stability (allowing applications with much 
larger time steps) may lead to a substantial improvement 
of the speed of MD simulations for magnetic liquids. It 
can also be used for lattices (then only Eqs. (4) and (5) 
must be employed) with arbitrary structures. These and 
related problems will be considered in a separate publi- 
cation. 
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FIG. 2. Transverse (a) and longitudinal (b) spin-spin, 
density-density (c), and spin-density (d) functions for a 
Heisenberg fluid versus frequency u>*. The curves correspond- 
ing to the wavenumbers k* = 1,2,3, and 5 are marked by " 1" , 
"2", "3", and "5", respectively. 

One peak can be identified for the function F^ s {k,ui) 
at each wavevector k. This peak is very sharp at small k 
and shifts to the right with increasing k. Such a quasipar- 
ticle behavior should be associated with the existence of 
transverse spin waves in the spin liquid. Up to three max- 
ima were observed for the component Fg S (k ) u>). While 
the first maximum at u) = corresponds to pure dif- 
fusive processes, the position of the second one coincides 
with that of the transverse spin wave peak, indicating the 
possibility of propagating longitudinal spin waves as well 
which, however, are damped much stronger. The origin 
of the third maximum in Fg S (k,u>) can be explained by 
the direct influence of the liquid subsystem on the spin 
one, because its position coincides with a peak position in 
F nn (k, bj). This last peak should be associated with prop- 
agative sound modes well established for liquid systems 
JL9[ , whereas a maximum of F nn (k,u>) at u> = repre- 
sents the well-known diffusive heat mode. The function 
F^ n (k,uj) behaves similarly to F nn (k,u}). In general the 
results obtained are in good agreement with the predic- 
tions of Ref . H . The additional possibility of longitudinal 
spin wave propagations in magnetic liquids at sound fre- 
quency can be considered as a new effect which has yet 
to be observed experimentally. Similar effect was found 
also (l8| in our MD calculations performed for model (1) 
at a higher temperature T > T c at the presence of an ex- 
ternal magnetic field. Taking into account the theoretical 
results of Ref. Q allows us to state that in the both cases 
considered the Brillouin sound peaks appear in F^ s (k,oj) 
due to magnetostriction caused by spin ordering. 
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